# pragma once

# include "cpp_tools_Headers.h"

namespace cpp_tools
{
    // 地球椭球体参数（WGS84）
    const double EARTH_A = 6378137.0;         // 长半轴，单位：米
    const double EARTH_F = 1.0 / 298.257223563; // 扁率
    const double EARTH_B = EARTH_A * (1 - EARTH_F); // 短半轴
    
    // 计算子午线弧长
    static double calculateMeridianArc(double B) {
        double e_sq = (EARTH_A*EARTH_A - EARTH_B*EARTH_B) / (EARTH_A*EARTH_A);
        double n = (EARTH_A - EARTH_B) / (EARTH_A + EARTH_B);
        
        double A = (EARTH_A + EARTH_B) / 2 * (1 + pow(n, 2)/4 + pow(n, 4)/64);
        double B_term = -3*n/2 + 9*pow(n, 3)/16 - 3*pow(n, 5)/32;
        double C_term = 15*pow(n, 2)/16 - 15*pow(n, 4)/32;
        double D_term = -35*pow(n, 3)/48 + 105*pow(n, 5)/256;
        double E_term = 315*pow(n, 4)/512;
        
        return A * (B + B_term*sin(2*B) + C_term*sin(4*B) + 
                    D_term*sin(6*B) + E_term*sin(8*B));
    }
    
    // 自动计算带号
    static int calculateZone(double longitude, bool is3Degree) {
        if (is3Degree) {
            // 3度带：从东经1.5°开始，每3°一带
            return static_cast<int>((longitude + 1.5) / 3.0) + 1;
        } else {
            // 6度带：从东经0°开始，每6°一带
            return static_cast<int>((longitude + 3.0) / 6.0) + 1;
        }
    }
    
    // 计算中央经线
    static double calculateCentralMeridian(int zone, bool is3Degree) {
        if (is3Degree) {
            return zone * 3.0;
        } else {
            return (zone - 1) * 6.0 + 3.0;
        }
    }
    
    /**
     * @brief 经纬度转高斯-克吕格坐标
     * @param latitude 纬度，单位：度
     * @param longitude 经度，单位：度
     * @param x 输出X坐标（米），包含带号
     * @param y 输出Y坐标（米），已加500km偏移
     * @param is3Degree 是否使用3度带（默认true）
     * @param zone 手动指定带号（默认0表示自动计算）
     */
    void latLonToGaussKruger(double latitude, double longitude, 
                            double& x, double& y, 
                            bool is3Degree = true, int zone = 0) {
        // 处理异常值
        if (latitude < -90 || latitude > 90 || longitude < -180 || longitude > 180) {
            std::cerr << "Invalid latitude or longitude!" << std::endl;
            x = y = 0;
            return;
        }
        
        // 自动计算带号（如果未指定）
        if (zone <= 0) {
            zone = calculateZone(longitude, is3Degree);
        }
        
        // 计算中央经线
        double L0 = calculateCentralMeridian(zone, is3Degree);
        
        // 转换为弧度
        double B = latitude * M_PI / 180.0;    // 纬度弧度
        double L = longitude * M_PI / 180.0;   // 经度弧度
        double L0_rad = L0 * M_PI / 180.0;     // 中央经线弧度
        
        // 计算经差
        double l = L - L0_rad;
        
        // 计算参数
        double e_sq = (EARTH_A*EARTH_A - EARTH_B*EARTH_B) / (EARTH_A*EARTH_A);
        double N = EARTH_A / sqrt(1 - e_sq * sin(B)*sin(B));  // 卯酉圈曲率半径
        double t = tan(B);
        double t_sq = t * t;
        double cosB = cos(B);
        double cosB_sq = cosB * cosB;
        double cosB_4 = cosB_sq * cosB_sq;
        double l_sq = l * l;
        double l_4 = l_sq * l_sq;
        
        // 计算子午线弧长
        double S = calculateMeridianArc(B);
        
        // 计算坐标
        x = S + N * t * cosB_sq * l_sq / 2 
              + N * t * (5 - t_sq + 9*e_sq*cosB_sq + 4*pow(e_sq, 2)*cosB_4) * cosB_4 * l_4 / 24;
        
        y = N * cosB * l 
              + N * (1 - t_sq + e_sq*cosB_sq) * cosB_sq * pow(l, 3) / 6 
              + N * (5 - 18*t_sq + pow(t, 4) + 14*e_sq*cosB_sq 
                     - 58*e_sq*t_sq*cosB_sq) * cosB_4 * pow(l, 5) / 120;
        
        // 加上带号和东偏值（500km）
        x = zone * 1000000 + x;  // 带号+坐标
        y += 500000;  // 东偏，避免Y出现负值
    }
    
    /**
     * @brief 高斯-克吕格坐标转经纬度
     * @param x 包含带号的X坐标（米）
     * @param y 已加500km偏移的Y坐标（米）
     * @param latitude 输出纬度，单位：度
     * @param longitude 输出经度，单位：度
     * @param is3Degree 是否使用3度带（默认true）
     */
    void gaussKrugerToLatLon(double x, double y, 
                            double& latitude, double& longitude, 
                            bool is3Degree = true) {
        // 提取带号
        int zone = static_cast<int>(x / 1000000);
        if (zone <= 0 || zone > (is3Degree ? 120 : 60)) {
            std::cerr << "Invalid zone number!" << std::endl;
            latitude = longitude = 0;
            return;
        }
        
        // 去除带号和东偏
        double X = x - zone * 1000000;
        y -= 500000;
        
        // 计算中央经线
        double L0 = calculateCentralMeridian(zone, is3Degree);
        
        // 计算底点纬度（迭代）
        double e_sq = (EARTH_A*EARTH_A - EARTH_B*EARTH_B) / (EARTH_A*EARTH_A);
        double Bf_rad = X / (EARTH_A * (1 - e_sq) * 
                         (1 + 3*e_sq/4 + 45*pow(e_sq, 2)/64 + 175*pow(e_sq, 3)/256));
        
        double delta;
        do {
            double sinBf = sin(Bf_rad);
            double term = 1 - e_sq * sinBf * sinBf;
            double M = EARTH_A * (1 - e_sq) / pow(term, 1.5);
            delta = (X - calculateMeridianArc(Bf_rad)) / M;
            Bf_rad += delta;
        } while (fabs(delta) > 1e-10);
        
        // 计算参数
        double t = tan(Bf_rad);
        double t_sq = t * t;
        double cosB = cos(Bf_rad);
        double cosB_sq = cosB * cosB;
        double cosB_4 = cosB_sq * cosB_sq;
        double N = EARTH_A / sqrt(1 - e_sq * sin(Bf_rad)*sin(Bf_rad));
        double y_sq = y * y;
        double y_4 = y_sq * y_sq;
        
        // 计算纬度
        double M = EARTH_A * (1 - e_sq) / pow(1 - e_sq * sin(Bf_rad)*sin(Bf_rad), 1.5);
        double B = Bf_rad - t/(2*M*N) * y_sq 
                  + t/(24*M*pow(N, 3)) * (5 + 3*t_sq + e_sq - 9*e_sq*t_sq) * y_4;
        
        // 计算经度
        double l = y/(N*cosB) - pow(y, 3)/(6*pow(N, 3)*cosB) * (1 + 2*t_sq + e_sq)
                  + pow(y, 5)/(120*pow(N, 5)*cosB) * (5 + 28*t_sq + 24*pow(t, 4) 
                  + 6*e_sq + 8*e_sq*t_sq);
        
        // 转换为度
        latitude = B * 180.0 / M_PI;
        longitude = (L0 * M_PI / 180.0 + l) * 180.0 / M_PI;
    }
} // namespace cpp_tools
